##------------------------------------------
##    policing the administrative state
##             dataverse code
##------------------------------------------


## replication code for the table and figures in the manuscript and appendix
## (also included is code for some of the key results in the appendix)


####   A. preamble

x <- read.csv(file = 'policing_dataverse.csv')
x <- subset(x, x$exec.agency == 1)

m <- aggregate(cbind(anprm, legal.dline,
                     interim.final, direct.final,
                     econ.significant, significant, sig, admin.time.trend,
                     ideology1, planned.first.year, planned.last.year, unified.gov,
                     ideology.gap) ~
                   agency.id + year.planned +
                   admin.nprm.planned,
               mean, data = x)
count <- aggregate(n.proposals ~
               agency.id + year.planned,
               sum, data = x)
m <- merge(m, count, by = c('agency.id', 'year.planned'))
a <- aggregate(audit ~ agency.id + year.planned,
               subset = nprm == 1,
               data = x, mean)
colnames(a)[3] <- 'inspection.rate'
m <- merge(m, a, by = c('agency.id', 'year.planned'))

# compliance
w <- aggregate(cbind(accept, returned, audit, nprm.stalled, nprm.planned) ~
                   year.planned + agency.id, sum,
               data = x, subset = audit == 1)
colnames(w)[colnames(w) == 'nprm.stalled'] <- c('nprm.stalled.audited')
colnames(w)[colnames(w) == 'nprm.planned'] <- c('nprm.planned.audited')
w <- merge(w, m, by = c('agency.id', 'year.planned'))
w$not.returned <- w$audit - w$returned
w$not.stalled <- w$nprm.planned.audited - w$nprm.stalled.audited

## withdrawals
x$withdrawn.before.nprm <- x$withdrawn
x$withdrawn.before.nprm[x$nprm == 1] <- 0
w2 <- aggregate(cbind(withdrawn.before.nprm, nprm.planned) ~
                   year.planned + agency.id, sum,
                data = x, subset = audit == 0)
colnames(w2)[colnames(w2) == 'nprm.planned'] <- c('nprm.planned.not.audited')
w <- merge(w, w2, by = c('agency.id', 'year.planned'))


####   B. main results


m1 <- glm(accept ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m2 <- glm(not.returned ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m3 <- glm(not.stalled ~ inspection.rate + offset(log(audit)) + #
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           ideology.gap + log(n.proposals) + sig,
                       data = w, family = quasipoisson())

m4 <- glm(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             ideology.gap + log(n.proposals) + sig,
                         data = w, family = quasipoisson())

# table for paper
k <- 4
n <- c('Inspection Rate', 'President-Agency Ideological Gap',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''

stargazer(m1, m2, m3, m4,
          type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          omit = 'agency|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'poisson1',
          title = 'Poisson Rate Model',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.year, blank))


####   C. robustness checks

####      C1. no controls


m1 <- glm(accept ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR'),
               data = w, family = quasipoisson())

m2 <- glm(not.returned ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR'),
               data = w, family = quasipoisson())

m3 <- glm(not.stalled ~ inspection.rate + offset(log(audit)) +
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR'),
                       data = w, family = quasipoisson())

m4 <- glm(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR'),
                         data = w, family = quasipoisson())

# table for paper
k <- 4
n <- c('Inspection Rate')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          style = 'apsr',
          omit = 'agency|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'poisson_no_control',
          title = 'Poisson Rate Model',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdraw Rate'),
          add.lines = list(blank, fe, fe.year, blank))



#####     C2. Time trends


m1 <- glm(accept ~ inspection.rate + offset(log(audit)) +
                   admin.time.trend +
                   admin.nprm.planned +
                   planned.first.year + planned.last.year +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + unified.gov + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m2 <- glm(not.returned ~ inspection.rate + offset(log(audit)) +
                   admin.time.trend +
                   admin.nprm.planned +
                   planned.first.year + planned.last.year +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + unified.gov + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m3 <- glm(not.stalled ~ inspection.rate + offset(log(audit)) +
                           admin.time.trend +
                           admin.nprm.planned +
                           planned.first.year + planned.last.year +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           ideology.gap + unified.gov + log(n.proposals) + sig,
                       data = w, family = quasipoisson())

m4 <- glm(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             admin.time.trend +
                             admin.nprm.planned +
                             planned.first.year + planned.last.year +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             ideology.gap + unified.gov + log(n.proposals) + sig,
                         data = w, family = quasipoisson())

# table for paper
k <- 4
n <- c('Inspection Rate', 'Administration Time Trend',
       'First Year of Administration', 'Last Year of Administration',
       'President-Agency Ideological Gap',
       'Unified Government', 'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          style = 'apsr',
          omit = 'agency|admin.np|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'timetrend',
          title = 'Administration Time Trend (Poisson Rate Model)',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.admin, blank))


#####     C3. Alternative agency ideology measures

## redo the aggregation from above using chen and johnson measures

## 1. chen johnson
## 2. intra-agency dispute

## redo preamble

m <- aggregate(cbind(sig, ideology1, chen.johnson, inter.agency.dispute,
                     ideology.gap) ~
               agency.id + year.planned,
               mean, data = x)
count <- aggregate(n.proposals ~
               agency.id + year.planned,
               sum, data = x)
m <- merge(m, count, by = c('agency.id', 'year.planned'))
a <- aggregate(audit ~ agency.id + year.planned,
               subset = nprm == 1,
               data = x, mean)
colnames(a)[3] <- 'inspection.rate'
m <- merge(m, a, by = c('agency.id', 'year.planned'))


w <- aggregate(cbind(accept, returned, audit, nprm.stalled, nprm.planned) ~
                   year.planned + agency.id, sum,
               data = x, subset = audit == 1)
colnames(w)[colnames(w) == 'nprm.stalled'] <- c('nprm.stalled.audited')
colnames(w)[colnames(w) == 'nprm.planned'] <- c('nprm.planned.audited')
w <- merge(w, m, by = c('agency.id', 'year.planned'))
w$not.returned <- w$audit - w$returned
w$not.stalled <- w$nprm.planned.audited - w$nprm.stalled.audited

x$withdrawn.before.nprm <- x$withdrawn
x$withdrawn.before.nprm[x$nprm == 1] <- 0
w2 <- aggregate(cbind(withdrawn.before.nprm, nprm.planned) ~
                   year.planned + agency.id, sum,
                data = x, subset = audit == 0)
colnames(w2)[colnames(w2) == 'nprm.planned'] <- c('nprm.planned.not.audited')
w <- merge(w, w2, by = c('agency.id', 'year.planned'))

## regressions

m1 <- glm(accept ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   chen.johnson + log(n.proposals) + sig,
               data = w, family = quasipoisson())


m2 <- glm(not.returned ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   chen.johnson + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m3 <- glm(not.stalled ~ inspection.rate + offset(log(audit)) +
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           chen.johnson + log(n.proposals) + sig,
                       data = w, family = quasipoisson())

m4 <- glm(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             chen.johnson + log(n.proposals) + sig,
                         data = w, family = quasipoisson())

## table for paper
k <- 4
n <- c('Inspection Rate', 'Ideological Gap (Chen-Johnson)',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          style = 'apsr',
          omit = 'agency|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'poisson1-chen',
          title = 'Chen-Johnson Ideology Measure (Poisson Rate Model)',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.year, blank))


## ------ inter-agency dispute


m1 <- glm(accept ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   inter.agency.dispute + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m2 <- glm(not.returned ~ inspection.rate + offset(log(audit)) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   inter.agency.dispute + log(n.proposals) + sig,
               data = w, family = quasipoisson())

m3 <- glm(not.stalled ~ inspection.rate + offset(log(audit)) +
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           inter.agency.dispute + log(n.proposals) + sig,
                       data = w, family = quasipoisson())

m4 <- glm(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             inter.agency.dispute + log(n.proposals) + sig,
                         data = w, family = quasipoisson())


k <- 4
n <- c('Inspection Rate', 'Inter-Agency Dispute',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          omit = 'agency.id|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'intra-agency',
          title = 'Intra-Agency Dispute (Poisson Rate Model)',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Abandon Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.year, blank))



#####     C4. linear models

## proposal-level summary stats

## get "x" from preamble

g <- summary(x[, c('interim.final', 'direct.final', 'legal.dline', 'econ.significant',
                   'significant', 'anprm')])
g <- g[-c(2,5,7),]
g[2,] <- g[4,]
g <- g[-4,]
g <- apply(g, 1, function(x) gsub('Min\\.|Max\\.|Mean|:', '', x))
g <- apply(g, 1, function(x) gsub(' +', '', x))
g <- apply(g, 1, function(x) substr(x, 1, 4))

n <- c('Interim Final Rule', 'Direct Final Rule', 'Legal Deadline Rule',
       'Economically Significant Rule', 'Significant Rule', 'Use of Advanced Notice Procedures')

rownames(g) <- n
sd1 <- sapply(x[, c('interim.final', 'direct.final', 'legal.dline', 'econ.significant',
                   'significant', 'anprm')], sd)
g <- cbind(g, round(sd1, 2))
xtable(g)

##  regression results

a <- aggregate(audit ~ agency.id + year.planned,
               subset = nprm == 1,
               data = x, mean)
colnames(a)[3] <- 'inspection.rate'
y <- merge(x, a, by = c('agency.id', 'year.planned'))

m1 <- lm(accept ~ inspection.rate +
             ideology.gap +
             relevel(as.factor(year.planned), ref = '2003') +
             relevel(as.factor(agency.id), ref = 'EPA-AR') +
             interim.final + direct.final + legal.dline +
             econ.significant + significant + anprm,
         data = y)

y$not.returned <- 1 - y$returned
m2 <- lm(not.returned ~ inspection.rate +
                  ideology.gap +
                  relevel(as.factor(year.planned), ref = '2003') +
                  relevel(as.factor(agency.id), ref = 'EPA-AR') +
                  interim.final + direct.final + legal.dline +
                  econ.significant + significant + anprm,
              data = y)

y$nprm.not.stalled <- 1 - y$nprm.stalled
m4 <- lm(nprm.not.stalled ~ inspection.rate +
             ideology.gap +
                  relevel(as.factor(year.planned), ref = '2003') +
                  relevel(as.factor(agency.id), ref = 'EPA-AR') +
                  interim.final + direct.final + legal.dline +
                  econ.significant + significant + anprm,
              data = y, subset = audit == 1)

y$withdrawn.before.nprm <- y$withdrawn
y$withdrawn.before.nprm[y$nprm == 1] <- 0
m5 <- lm(withdrawn ~ inspection.rate +
                  ideology.gap +
                  relevel(as.factor(year.planned), ref = '2003') +
                  relevel(as.factor(agency.id), ref = 'EPA-AR') +
                  interim.final + direct.final + legal.dline +
                  econ.significant + significant + anprm,
              data = y, subset = audit == 0)


# clustered standard errors
cluster <- function(dat, fm, cluster){
           require(sandwich, quietly = TRUE)
           require(lmtest, quietly = TRUE)
           M <- length(unique(cluster))
           N <- length(cluster)
           K <- fm$rank
           dfc <- (M/(M-1))*((N-1)/(N-K))
           uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
           uj <- na.omit(uj)
           vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
           coeftest(fm, vcovCL) }


y1 <- data.frame(with(y, na.omit(cbind(accept, not.returned,
                                       inspection.rate,
                                       ideology.gap,
                                       interim.final, direct.final, legal.dline,
                                       econ.significant,
                                       significant, anprm,
                                       agency.id,
                                       year.planned))))
y1$cluster <- paste(y1$agency.id, y1$year.planned, sep = '-')
se1 <- cluster(y1, m1, y1$cluster)
se2 <- cluster(y1, m2, y1$cluster)
se4 <- cluster(y1, m4, y1$cluster)

y2 <- data.frame(with(y, na.omit(cbind(nprm.stalled, inspection.rate, admin.nprm.planned,
                                       ideology.gap,
                                       interim.final, direct.final, legal.dline,
                                       econ.significant,
                                       significant, anprm,
                                       agency.id,
                                       year.planned, audit))))
y2 <- subset(y2, audit == 0)
y2$cluster <- paste(y2$agency.id, y2$year.planned, sep = '-')

se5 <- cluster(y2, m5, y2$cluster)


# table for paper
k <- 6
n <- c('Inspection Rate', 'President-Agency Ideological Gap',
       'Interim Final Rule', 'Direct Final Rule', 'Legal Deadline Rule',
       'Economically Significant Rule', 'Significant Rule', 'Use of Advanced Notice Procedures')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Administration Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
controls <- c('Rule-Type Controls', rep(c('', '\\checkmark'), 3))
weights <- c('Weights for \\# of Proposals', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m4, m5,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          style = 'apsr',
          omit = 'agency|year|admin|Constant',
          covariate.labels = n,
          no.space = T,
          se = list(se1[,2], se2[,2], se4[,2], se5[,2]),
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'lpm',
          title = 'Linear Probability Model',
          dep.var.labels = c('Change Rate', 'Return Rate',
            'Withdrawn Rate', 'Stall Rate', 'Withdraw Rate'),
          add.lines = list(blank, fe, fe.year, blank))



##  --- log linear model


## uses w matrix from 'main analysis'

m1 <- lm(log(accept+1) ~ inspection.rate + log(audit+1) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + log(n.proposals) + sig,
               data = w)

m2 <- lm(log(not.returned+1) ~ inspection.rate + log(audit) +
                   relevel(as.factor(year.planned), ref = '2003') +
                   relevel(as.factor(agency.id), ref = 'EPA-AR') +
                   ideology.gap + log(n.proposals) + sig,
               data = w)

m3 <- lm(log(not.stalled+1) ~ inspection.rate + log(audit) +
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           ideology.gap + log(n.proposals) + sig,
                       data = w)

m4 <- lm(log(withdrawn.before.nprm+1) ~ inspection.rate + log(nprm.planned.not.audited) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             ideology.gap + log(n.proposals) + sig,
                         data = w)

# table for paper
k <- 4
n <- c('Inspection Rate', 'N', 'N2', 'President-Agency Ideological Gap',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''


stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          style = 'apsr',
          omit = 'agency|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'ols2',
          title = 'Main Results Using Log-Linear Model',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.year, blank))


#####     B1. GEE Poisson

library(gee)

w$agency.id <- droplevels(w$agency.id)

## exchangeable correlation structure: cor(y_it, y_it') = a

m1 <- gee(accept ~ inspection.rate + offset(log(audit)) +
              relevel(as.factor(year.planned), ref = '2003') +
              relevel(as.factor(agency.id), ref = 'EPA-AR') +
              ideology.gap + log(n.proposals) + sig,
          data = w, family = 'poisson', id = year.planned,
          corstr = 'exchangeable', maxiter = 1000)

m2 <- gee(not.returned ~ inspection.rate + offset(log(audit)) +
              relevel(as.factor(year.planned), ref = '2003') +
              relevel(as.factor(agency.id), ref = 'EPA-AR') +
              ideology.gap + log(n.proposals) + sig,
          data = w, family = 'poisson', id = year.planned,
          corstr = 'exchangeable', maxiter = 1000)

m3 <- gee(not.stalled ~ inspection.rate + offset(log(audit)) +
                           relevel(as.factor(year.planned), ref = '2003') +
                           relevel(as.factor(agency.id), ref = 'EPA-AR') +
                           ideology.gap + log(n.proposals) + sig,
          data = w, family = 'poisson', id = year.planned,
          corstr = 'exchangeable', maxiter = 1000)

m4 <- gee(withdrawn.before.nprm ~ inspection.rate + offset(log(nprm.planned.not.audited)) +
                             relevel(as.factor(year.planned), ref = '2003') +
                             relevel(as.factor(agency.id), ref = 'EPA-AR') +
                             ideology.gap + log(n.proposals) + sig,
          data = w, family = 'poisson', id = year.planned,
          corstr = 'exchangeable', maxiter = 1000)

# table for paper
k <- 4
n <- c('Inspection Rate', 'President-Agency Ideological Gap',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe <- c('Agency Fixed Effects', rep(c('\\checkmark'), k))
fe.admin <- c('Admin Fixed Effects', rep(c('\\checkmark'), k))
fe.year <- c('Year Fixed Effects', rep(c('\\checkmark'), k))
blank <- ''

stargazer(m1, m2, m3, m4,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          omit = 'agency|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'poisson1',
          title = 'Poisson Rate Model',
          dep.var.labels = c('Accept Rate', 'Return Rate',
            'Stall Rate', 'Withdrawn Rate'),
          add.lines = list(blank, fe, fe.year, blank))



#####     C2. non-parametric estimation


audit1 <- lm(inspection.rate ~
                 as.factor(year.planned) + log(n.proposals) +
                 ideology.gap +
                 sig + agency.id,
                 data = w, weights = audit)


compliance1 <- lm(I(accept/audit) ~ log(n.proposals) +
                 as.factor(year.planned) +
                     ideology.gap +
                     sig +
                 agency.id,
              data = w, weights = audit)

compliance2 <- lm(I(not.returned/audit) ~ log(n.proposals) +
                 as.factor(year.planned) +
                     ideology.gap +
                     sig +
                 agency.id,
              data = w, weights = audit)

compliance3 <- lm(I(not.stalled/nprm.planned.audited) ~ log(n.proposals) +
                 as.factor(year.planned) +
                     ideology.gap +
                     sig +
                 agency.id,
              data = w, weights = audit)

compliance4 <- lm(I(withdrawn.before.nprm/nprm.planned.not.audited) ~ log(n.proposals) +
                 as.factor(year.planned) +
                     ideology.gap +
                     sig +
                 agency.id,
              data = w, weights = audit)


cor.test(compliance1$resid, audit1$resid, method = 'kendall')
cor.test(compliance2$resid, audit1$resid, method = 'kendall')
cor.test(compliance3$resid, audit1$resid, method = 'kendall')
cor.test(compliance4$resid, audit1$resid, method = 'kendall')



#####     C3. sensitivity analysis

## run preamble to 'main analysis' first

## test 1, accept rate
d <- c()
for(i in 1:1000) {
    n <- sample(1:nrow(w), .7*nrow(w), replace = F)
    v <- w[n,]
    m <- glm(accept ~ inspection.rate + offset(log(audit)) +
             as.factor(year.planned) + ideology.gap +
             log(n.proposals) + sig + agency.id,
             data = v, family = quasipoisson())
    a <- data.frame(est = summary(m)$coef['inspection.rate',1],
                    se = summary(m)$coef['inspection.rate',2])
    d <- rbind(d, a)
}
mean(as.numeric(abs(d$est) > 1.96*d$se)) # .26
mean(as.numeric(d$est < 0)) # 0

## test 2, return rate
d <- c()
for(i in 1:1000) {
    n <- sample(1:nrow(w), .7*nrow(w), replace = F)
    v <- w[n,]
    m <- glm(not.returned ~ inspection.rate + offset(log(1+audit)) +
             as.factor(year.planned) + ideology.gap +
             log(n.proposals) + sig +
             agency.id, data = v, family = quasipoisson())
    a <- data.frame(est = summary(m)$coef['inspection.rate',1],
                    se = summary(m)$coef['inspection.rate',2])
    d <- rbind(d, a)
}
mean(as.numeric(abs(d$est) > 1.96*d$se)) # 1


## test 3, stall rate
d <- c()
for(i in 1:1000) {
    n <- sample(1:nrow(w), .7*nrow(w), replace = F)
    v <- w[n,]
    m <- glm(not.stalled ~ offset(log(audit)) + inspection.rate +
             as.factor(year.planned) + ideology.gap +
             log(n.proposals) + sig +
             agency.id, data = v, family = quasipoisson())
    a <- data.frame(est = summary(m)$coef['inspection.rate',1],
                    se = summary(m)$coef['inspection.rate',2])
    d <- rbind(d, a)
}
mean(as.numeric(abs(d$est) > 1.96*d$se)) # 1


## test 4, withdrawal rate
d <- c()
for(i in 1:1000) {
    n <- sample(1:nrow(w), .7*nrow(w), replace = F)
    v <- w[n,]
    m <- glm(withdrawn.before.nprm ~ offset(log(nprm.planned.not.audited)) + inspection.rate +
             as.factor(year.planned) + ideology.gap +
             log(n.proposals) + sig + agency.id,
             data = v, family = quasipoisson())
    a <- data.frame(est = summary(m)$coef['inspection.rate',1],
                    se = summary(m)$coef['inspection.rate',2])
    d <- rbind(d, a)
}
mean(as.numeric(abs(d$est) > 1.96*d$se)) # 1


#####     E1. executive branch vs. IRC


x <- read.csv(file = '~/Dropbox/research/deterrence/jop_files/policing_dataverse.csv')

## aggregate to agency-year level
m <- aggregate(cbind(anprm, legal.dline,
                     interim.final, direct.final,
                     econ.significant, significant,sig,
                     ideology.gap, planned.first.year, planned.last.year,
                     exec.agency, unified.gov) ~
                   agency.id + year.planned + admin.nprm.planned,
               mean,
               data = x)

count <- aggregate(n.proposals ~
               agency.id + year.planned,
               sum, data = x)
m <- merge(m, count, by = c('agency.id', 'year.planned'))

x$withdrawn.before.nprm <- x$withdrawn
x$withdrawn.before.nprm[x$nprm == 1] <- 0
w2 <- aggregate(cbind(withdrawn.before.nprm, nprm.planned) ~
                   year.planned + agency.id, sum,
                data = x, subset = audit == 0)
colnames(w2)[colnames(w2) == 'nprm.planned'] <- c('nprm.planned.not.audited')
w2 <- merge(w2, m, by = c('agency.id', 'year.planned'))


## regresssions

library(gee)
w$agency.id <- droplevels(w$agency.id)

m2 <- gee(withdrawn.before.nprm ~ exec.agency +
               offset(log(nprm.planned.not.audited)),
                  data = w2,
          family = 'poisson',
          id = agency.id, corstr = 'exchangeable', maxiter = 1000)

m2.b <- gee(withdrawn.before.nprm ~
               offset(log(nprm.planned.not.audited)) +
                  exec.agency + year.planned +
               ideology.gap +
               log(n.proposals) + sig,
            data = w2, family = 'poisson',
            id = agency.id, corstr = 'exchangeable', maxiter = 1000)

## matching analysis
library(MatchIt)

match2 <- matchit(exec.agency ~  sig + ideology.gap + n.proposals + year.planned,
                  data = w2, method = 'nearest')
z2 <- match.data(match2)
m2.matched <- update(m2.b, data = z2)


# table for paper
k <- 3
n <- c('Executive Branch Agency', 'President-Agency Ideological Gap',
       'N. of Proposals (logged)',
       'Prop. of Proposals Significant')
fe.year <- c('Year Fixed Effects', '', '\\checkmark', '\\checkmark')
blank <- ''


stargazer(m2, m2.b, m2.matched,
          ## type = 'text',
          omit.stat=c("f", 'ser', 'adj.rsq'),
          omit = 'admin|year|Constant',
          covariate.labels = n,
          no.space = T,
          star.cutoffs = c(.05, .01, .001),
          font.size = 'footnotesize',
          label = 'irc',
          title = 'Comparison with Indepenent Commissions (Poisson Rate Model)',
          add.lines = list(blank, fe.year, blank))


####   D. figures & tables


####      D1. summary statistics

## (run after running the preamble to the Poisson models above)

library(xtable)

length(unique(x$year.planned))  # 20 years
length(unique(x$agency.id))  # 42 agencies

## 1. summary of variables

# create rates
w$compliance.rate1 <- w$accept / w$audit
w$compliance.rate2 <- w$not.returned / w$audit
w$compliance.rate3 <- w$not.stalled / w$audit
w$withdraw.rate <- w$withdrawn.before.nprm / w$nprm.planned.not.audited

y <- summary(w[, c('inspection.rate', 'compliance.rate1', 'compliance.rate2',
                   'compliance.rate3', 'withdraw.rate',
                   'year.planned', 'sig', 'ideology.gap', 'n.proposals')])
y <- y[-c(2,5,7),]
y[2,] <- y[4,]
y <- y[-4,]
y <- apply(y, 1, function(x) gsub('Min\\.|Max\\.|Mean|:', '', x))
y <- apply(y, 1, function(x) gsub(' +', '', x))
y <- apply(y, 1, function(x) substr(x, 1, 4))

n <- c('Inspection Rate', 'Compliance Rate 1', 'Compliance Rate 2', 'Compliance Rate 3',
       'Withdrawal Rate', 'Year', 'Prop. Significant',
       'Ideology Gap', 'Number of Proposals')

rownames(y) <- n
sd1 <- sapply(w[, c('inspection.rate', 'compliance.rate1', 'compliance.rate2',
                    'compliance.rate3', 'withdraw.rate',
                     'year.planned', 'sig', 'ideology.gap', 'n.proposals')], sd)
y <- cbind(y, round(sd1, 2))
xtable(y)



#####     D2. figure: EPA audit and hit rate plots



w$accept.rate <- w$accept / w$audit
w$hit.rate <- 1 - w$accept.rate
w$agency.admin <- paste(w$admin.nprm.planned, w$agency.id, sep = '-')


v <- data.frame(rate = c(w$inspection.rate, w$hit.rate),
                Rates = c(rep('Inspection Rate', nrow(w)), rep('Hit Rate', nrow(w)) ),
                year = rep(w$year.planned, 2),
                agency = rep(w$agency.admin, 2))

v$Rates <- factor(v$Rates, levels = c('Inspection Rate','Hit Rate'))

v <- subset(v, agency == 'clinton-EPA-AR')

v1 <- v[v$Rates == 'Inspection Rate',]
v2 <- v[v$Rates == 'Hit Rate',]

par(mfrow = c(1,1), mar = c(3, 3, 2, 1))
plot(v$year, v$rate,
     xlab = '', ylim = c(0, 1), cex.axis = 1.4, las = 1,
     ylab = '',
     main = '')
points(v1$year, v2$rate, pch = 19)
abline(lm(rate ~ year, data = v1))
abline(lm(rate ~ year, data = v2))
legend('bottomleft', legend = c('Hit Rate', 'Inspection Rate'), pch = c(19, 1), cex = 1.3)


#####     D3. regression result plot (year fixed effects)

## run 'main results' first (section B, models: m1, m2, m3, m4)

## note: the simulated confidence regions say nothing about the
##       statistical significance of the relationship; they only give
##       a sense of the uncertainty in the trend line

library(arm)

## pdf('~/Dropbox/research/deterrence/jop_files/inspection_curves_year_effects.pdf',
##     height = 4)

par(mfrow = c(1,2), mar = c(5, 4, 2, 1))


# panel 1
fit2 <- sim(m2)
plot(0:1, c(.5,1), col = 'white',
     xlab = 'Inspection Rate',
     ylab = 'Compliance Rate',
     main = 'Deterrence Test 1')
for(i in 1:30) curve( exp(coef(fit2)[i,'(Intercept)'] +
                           coef(fit2)[i,'inspection.rate'] *x) ,
                     add = TRUE, col = 'lightgrey')
curve( exp(m2$coef['(Intercept)'] +
           m2$coef['inspection.rate'] *x) ,
      add = TRUE)
rug(jitter(w$inspection.rate, amount = .03), lwd = .2)

# panel 2
fit5 <- sim(m4)
plot(0:1, c(0,.5), col = 'white',
     xlab = 'Inspection Rate',
     ylab = 'Withdrawal Rate',
     main = 'Deterrence Test 2')
for(i in 1:20) curve( exp(coef(fit5)[i, '(Intercept)'] +
                           coef(fit5)[i,'inspection.rate'] *x) ,
                     add = TRUE, col = 'lightgrey')
curve( exp(m4$coef['(Intercept)'] +
           m4$coef['inspection.rate'] *x) ,
      add = TRUE)
rug(jitter(w$inspection.rate, amount = .03), lwd = .2)

## dev.off()


#####     D4. table of agencies

## run the preamble first (part A)

w <- w[!duplicated(w$agency.id),]

a <- subset(x, select = c(agency.id, agency, bureau) )
a <- a[!duplicated(a),]
y <- merge(w, a, by = c('agency.id'), all.x = T)

y <- aggregate(cbind(inspection.rate, n.proposals) ~ agency + agency.id + bureau, data = y, mean)

y$agency <- as.character(y$agency)
y$bureau <- as.character(y$bureau)

dept <- gsub('(.*)-.*', '\\1', y$agency.id)
agency <- gsub('.*-(.*)', '\\1', y$agency.id)
y$agency.id <- as.character(y$agency.id)
y$agency.id[dept == agency] <- dept[dept == agency]

y <- y[,c(1,3:5)]
colnames(y) <- c('Department','Agency','Audit Rate', 'N. Proposed')

y$Department[y$Department == y$Agency] <- 'x'
y <- y[order(y$Department),]
y$Department[y$Department == 'x'] <- ''
rownames(y) <- 1:nrow(y)
mean(y[,3])
sum(y[,4])

print.xtable(xtable(y))
